03. CTmax_analysis_resp_analysis
Elliott Schmidt
2024-08-01
Scenario
Within this chapter developmental plasticity was explored within Acanthochromis polyacanthus that were collected from two different regions (i.e., low-latitude, Cairns, and high-latitude, Mackay). Fish were held in common garden experiments at 28.5 C. Clutch data was collected along with parental morphometeric data to determine how fish from each population performed at 28.5 C, a temperature that was shown to produce similar absoluate aerobic scope performances in a previous study [link]. Once hatched offspring were placed into three different temperature treatments, 28.5, 30, and 31.5 C. After approximately 5-6 months offspring length and weight was measured, as well as CTmax and respiration during CTmax trials.
Load packages
library(tidyverse) # data manipulation
library(ggpubr) # figure arrangement
library(brms) # Bayesian models
library(StanHeaders)# needed to run Bayesian models
library(rstan) # needed to run Bayesian models
library(standist) # needs to be installed
library(bayesplot) # needed for MCMC diagnostics
library(DHARMa) # model validation
library(ggdist) # partial plots
library(tidybayes) # partial plots
library(broom.mixed) # model investigation
library(emmeans) # pairwise comparisons
library(rstanarm) # pairwise comparisons - need for emmeans
library(ggpubr) # arranging figures
library(mgcv) # GAM
library(modelr) # plotting
library(segmented) # piece wise regression
library(strucchange) #piece wise regressionImport data
This data set has no point outliers however, in the next chunk 3 samples will be discarded due to poor data quailty
Data manipulation
rows_of_data <- count(lat_resp_dat, sampleID)
lat_resp_dat2 <- lat_resp_dat |>
mutate(dev.temp = as.factor(dev.temp),
replicate = str_sub(sampleID, -1,-1),
population = factor(population)) |>
# number of observations = 5758
filter(sampleID != "56.CARL.137.28,5.1", # 5777 - 76 = 5682
sampleID != "56.CARL.137.28,5.2", # 5701 - 64 = 5618
sampleID != "60.LCKM.152.30.1" # 5637 - 76 = 5542
) |>
filter(time_lag_sec >2001) |> # remove samples from first 5 cycles (i.e., first 33 minutes)
group_by(sampleID) |>
mutate(max_value_index = which.max(rate.output),
row_number = row_number(),
max.rate = max(rate.output),
low.rate = mean(rate.output[order(rate.output)[1:10]]),
net.rate = max.rate - low.rate) |>
filter(row_number <= max_value_index) |>
ungroup() |>
mutate(region = factor(region),
dev.temp =factor(dev.temp),
level = as.factor(case_when(tank >= 1 & tank <= 199 ~ 1,
tank >= 200 & tank <= 299 ~ 2,
tank >= 300 & tank <= 399 ~ 3,
TRUE ~ NA_real_)),
female = factor(female)) |>
unite("dev.temp.region",c(dev.temp,region), remove=FALSE) |>
mutate(dev.temp.region = factor(dev.temp.region))Exploratory data analysis
2nd order polynomial
eda1 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output)) +
geom_point(alpha=0.5, color = "black") +
geom_smooth(method="lm",
se=TRUE,
fill=NA,
formula = y~poly(x, 2, raw=TRUE), color = "red") +
theme_classic() +
ggtitle("All data points - 2nd order polynomial")
eda2 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=region)) +
geom_point(alpha=0.5) +
geom_smooth(method="lm",
se=TRUE,
fill=NA,
formula = y~poly(x, 2, raw=TRUE), color = "red") +
facet_wrap(~ region) + theme_classic() +
theme(legend.position = "none") +
ggtitle("Region - 2nd order polynomial")
eda3 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=dev.temp)) +
geom_point(alpha=0.1) +
geom_smooth(method="lm",
se=TRUE,
fill=NA,
formula = y~poly(x, 2, raw=TRUE)) +
facet_wrap(~ region) + theme_classic() +
theme(legend.position = "none") +
ggtitle("Region - 2nd order polynomial")
eda4 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=region)) +
geom_point(alpha=0.1) +
geom_smooth(method="lm",
se=TRUE,
fill=NA,
formula = y~poly(x, 2, raw=TRUE)) +
facet_wrap(~ dev.temp) + theme_classic() +
theme(legend.position = "none") +
ggtitle("Region - 2nd order polynomial")
eda5 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=population)) +
geom_point(alpha=0.5) +
geom_smooth(method="lm",
se=TRUE,
fill=NA,
formula = y~poly(x, 2, raw=TRUE), color = "red") +
facet_wrap(~ population) + theme_classic() +
theme(legend.position = "none") +
ggtitle("Population - 2nd order polynomial")
eda6 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=dev.temp)) +
geom_point(alpha=0.1) +
geom_smooth(method="lm",
se=TRUE,
fill=NA,
formula = y~poly(x, 2, raw=TRUE)) +
facet_wrap(~ population) + theme_classic() +
theme(legend.position = "none") +
ggtitle("Population - 2nd order polynomial")
eda.fig <- ggarrange(eda1,eda2,eda3,eda4,eda5,eda6,
nrow = 6,
ncol=1); eda.fig
3rd order polynomial
eda1 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output)) +
geom_point(alpha=0.5, color = "black") +
geom_smooth(method="lm",
se=TRUE,
fill=NA,
formula = y~poly(x, 3, raw=TRUE), color = "red") +
theme_classic() +
ggtitle("All data points - 2nd order polynomial")
eda2 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=region)) +
geom_point(alpha=0.5) +
geom_smooth(method="lm",
se=TRUE,
fill=NA,
formula = y~poly(x, 3, raw=TRUE), color = "red") +
facet_wrap(~ region) + theme_classic() +
theme(legend.position = "none") +
ggtitle("Region - 2nd order polynomial")
eda3 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=dev.temp)) +
geom_point(alpha=0.1) +
geom_smooth(method="lm",
se=TRUE,
fill=NA,
formula = y~poly(x, 3, raw=TRUE)) +
facet_wrap(~ region) + theme_classic() +
theme(legend.position = "none") +
ggtitle("Region - 2nd order polynomial")
eda4 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=region)) +
geom_point(alpha=0.1) +
geom_smooth(method="lm",
se=TRUE,
fill=NA,
formula = y~poly(x, 3, raw=TRUE)) +
facet_wrap(~ dev.temp) + theme_classic() +
theme(legend.position = "none") +
ggtitle("Region - 2nd order polynomial")
eda5 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=population)) +
geom_point(alpha=0.5) +
geom_smooth(method="lm",
se=TRUE,
fill=NA,
formula = y~poly(x, 3, raw=TRUE), color = "red") +
facet_wrap(~ population) + theme_classic() +
theme(legend.position = "none") +
ggtitle("Population - 2nd order polynomial")
eda6 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=dev.temp)) +
geom_point(alpha=0.1) +
geom_smooth(method="lm",
se=TRUE,
fill=NA,
formula = y~poly(x, 3, raw=TRUE)) +
facet_wrap(~ population) + theme_classic() +
theme(legend.position = "none") +
ggtitle("Population - 2nd order polynomial")
eda.fig <- ggarrange(eda1,eda2,eda3,eda4,eda5,eda6,
nrow = 6,
ncol=1); eda.fig
Fit model [random factors]
First we will place random factors only within out model and see if they do a good job explaining the variance within our model. We will also be include priors which will be based off of out length data.
Hypothesis test will include:
- Null model
- EGG_COUNT ~ 1 + (1| FEMALE) + (1| TANK)
- EGG_COUNT ~ 1 + (1| FEMALE) + (1| TANK) + (1| LEVEL)
- EGG_COUNT ~ 1 + (1| FEMALE) + (1| TANK) + (1| LEVEL) + (1| POPULATION)
- EGG_COUNT ~ 1 + (1| FEMALE) + (1| TANK) + (1| LEVEL) + (1| POPULATION)+ (1| CLUTCH_ORDER)
lat_resp_dat2 |>
group_by(region,dev.temp) |>
summarise(mean(na.omit(rate.output)),
sd(na.omit(rate.output))) ## `summarise()` has grouped output by 'region'. You can override using the
## `.groups` argument.Priors
Models
Null
f.model.null <- bf(rate.output2 ~ 1,
family = gaussian())
model.null <- brm(f.model.null,
data = lat_resp_dat2,
prior = rate.priors2,
warmup = 500,
iter = 5000,
seed=123,
cores=2,
save_pars = save_pars(all=TRUE),
sample_prior = "yes",
chains = 2,
thin = 5,
control = list(adapt_delta=0.95)) ## Compiling Stan program...## Start samplingModel1
f.model1 <- bf(rate.output2 ~ 1 + (1| female) + (1| tank),
family=gaussian())
model1 <- brm(f.model1,
data = lat_resp_dat2,
prior = rate.priors2,
warmup = 500,
iter = 5000,
seed=123,
cores=2,
save_pars = save_pars(all=TRUE),
sample_prior = "yes",
chains = 2,
thin = 5,
control = list(adapt_delta=0.95)) ## Compiling Stan program...## Start samplingModel2
f.model2 <- bf(rate.output2 ~ 1 + (1| female) + (1| tank) + (1| level),
family=gaussian())
model2 <- brm(f.model2,
data = lat_resp_dat2,
prior = rate.priors2,
warmup = 500,
iter = 5000,
seed=123,
cores=2,
save_pars = save_pars(all=TRUE),
sample_prior = "yes",
chains = 2,
thin = 5,
control = list(adapt_delta=0.95)) ## Compiling Stan program...## Start sampling## Warning: There were 26 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.## Warning: There were 753 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
## https://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded## Warning: Examine the pairs() plot to diagnose sampling problemsModel3
f.model3 <- bf(rate.output2 ~ 1 + (1| female) + (1| tank) + (1| level) + (1| population),
family=gaussian())
model3 <- brm(f.model3,
data = lat_resp_dat2,
prior = rate.priors2,
warmup = 500,
iter = 5000,
seed=123,
cores=2,
save_pars = save_pars(all=TRUE),
sample_prior = "yes",
chains = 2,
thin = 5,
control = list(adapt_delta=0.95)) ## Compiling Stan program...## Start sampling## Warning: There were 18 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.## Warning: There were 85 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
## https://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded## Warning: Examine the pairs() plot to diagnose sampling problemsModel4
f.model4 <- bf(rate.output2 ~ 1 + (1| female) + (1 | tank) + (1| level) + (1| population)+ (1| clutch_order),
family=gaussian())
model4 <- brm(f.model4,
data = lat_resp_dat2,
prior = rate.priors2,
warmup = 500,
iter = 5000,
seed=123,
cores=2,
save_pars = save_pars(all=TRUE),
sample_prior = "yes",
chains = 2,
thin = 5,
control = list(adapt_delta=0.95)) ## Compiling Stan program...## Start sampling## Warning: There were 74 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.## Warning: There were 14 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
## https://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded## Warning: Examine the pairs() plot to diagnose sampling problemsLOO
## Output of model 'model.null':
##
## Computed from 1800 by 4776 log-likelihood matrix.
##
## Estimate SE
## elpd_loo -1446.4 81.9
## p_loo 3.9 0.3
## looic 2892.8 163.8
## ------
## MCSE of elpd_loo is 0.0.
## MCSE and ESS estimates assume MCMC draws (r_eff in [1.0, 1.0]).
##
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
##
## Output of model 'model1':
##
## Computed from 1800 by 4776 log-likelihood matrix.
##
## Estimate SE
## elpd_loo 91.4 83.1
## p_loo 56.1 2.3
## looic -182.8 166.2
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.8, 1.1]).
##
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
##
## Output of model 'model2':
##
## Computed from 1800 by 4776 log-likelihood matrix.
##
## Estimate SE
## elpd_loo 91.3 83.2
## p_loo 56.3 2.3
## looic -182.7 166.4
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.8, 1.1]).
##
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
##
## Output of model 'model3':
##
## Computed from 1800 by 4776 log-likelihood matrix.
##
## Estimate SE
## elpd_loo 91.3 83.1
## p_loo 56.3 2.3
## looic -182.5 166.3
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.6, 1.1]).
##
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
##
## Output of model 'model4':
##
## Computed from 1800 by 4776 log-likelihood matrix.
##
## Estimate SE
## elpd_loo 90.8 83.2
## p_loo 56.8 2.3
## looic -181.7 166.4
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.3, 1.1]).
##
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
##
## Model comparisons:
## elpd_diff se_diff
## model1 0.0 0.0
## model2 -0.1 0.4
## model3 -0.1 0.4
## model4 -0.5 0.5
## model.null -1537.8 77.4Fit model [fixed factors - polynomials]
Prior investigation
lat_resp_dat2 |>
group_by(region,dev.temp) |>
summarise(mean(na.omit(rate.output2)),
sd(na.omit(rate.output2))) ; lat_resp_dat2 |>
group_by(region,dev.temp) |>
summarise(mean(na.omit(rate.output)),
sd(na.omit(rate.output))) ## `summarise()` has grouped output by 'region'. You can override using the
## `.groups` argument.## `summarise()` has grouped output by 'region'. You can override using the
## `.groups` argument.Priors
# change to reflect rate.output2
rate.priors2 <- prior(normal(0.635, 0.20), class="Intercept") +
prior(normal(0, 0.20), class="b") +
prior(student_t(3, 0, 0.45), class = "sigma")
rate.priors <- prior(normal(430, 0.25), class="Intercept") +
prior(normal(0, 40), class="b") +
prior(student_t(3, 0, 195), class = "sigma")Models
Model1.p1
Model1.p1
f.model1.p1 <- bf(rate.output2 ~ 1 +
dev.temp*region + Temperature +
scale(mass, center=TRUE, scale=TRUE) +
(1| female) + (1| tank) ,
family=gaussian())
model1.p1 <- brm(f.model1.p1,
data = lat_resp_dat2,
prior = rate.priors2,
warmup = 500,
iter = 5000,
seed=123,
cores=2,
save_pars = save_pars(all=TRUE),
sample_prior = "yes",
chains = 2,
thin = 5,
control = list(adapt_delta=0.95)) ## Compiling Stan program...## Start samplingDHARMa
#--- distribution check ---#
pp_check(model1.p1, type = 'dens_overlay_grouped', ndraws=150, group="region")
#--- DHARMa checks ---#
preds <- posterior_predict(model1.p1, ndraws=250, summary=FALSE)
model1.p1_resids <- createDHARMa(simulatedResponse = t(preds),
observedResponse = lat_resp_dat2$rate.output,
fittedPredictedResponse = apply(preds, 2, mean),
integerResponse = 'student')
plot(model1.p1_resids) ; testDispersion(model1.p1_resids)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 972987, p-value < 2.2e-16
## alternative hypothesis: two.sidedMCMC diagnostics
## 'pars' not specified. Showing first 10 parameters by default.
## 'pars' not specified. Showing first 10 parameters by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## 'pars' not specified. Showing first 10 parameters by default.
model1.p2
model1.p2
f.model1.p2 <- bf(rate.output2 ~ 1 +
dev.temp*region + poly(Temperature, 2) +
scale(mass, center=TRUE, scale=TRUE) +
(1| female) + (1| tank) ,
family=gaussian())
model1.p2 <- brm(f.model1.p2,
data = lat_resp_dat2,
prior = rate.priors2,
warmup = 500,
iter = 5000,
seed=123,
cores=2,
save_pars = save_pars(all=TRUE),
sample_prior = "yes",
chains = 2,
thin = 5,
control = list(adapt_delta=0.95)) ## Compiling Stan program...## Start samplingDHARMa
#--- distribution check ---#
pp_check(model1.p2, type = 'dens_overlay_grouped', ndraws=150, group="region")
#--- DHARMa checks ---#
preds <- posterior_predict(model1.p2, ndraws=250, summary=FALSE)
model1.p2_resids <- createDHARMa(simulatedResponse = t(preds),
observedResponse = lat_resp_dat2$rate.output,
fittedPredictedResponse = apply(preds, 2, mean),
integerResponse = 'student')
plot(model1.p2_resids) ; testDispersion(model1.p2_resids)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1044481, p-value < 2.2e-16
## alternative hypothesis: two.sidedMCMC diagnostics
## 'pars' not specified. Showing first 10 parameters by default.
## 'pars' not specified. Showing first 10 parameters by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## 'pars' not specified. Showing first 10 parameters by default.
model1.p3
model1.p3
f.model1.p3 <- bf(rate.output2 ~ 1 +
dev.temp*region + poly(Temperature, 3) +
scale(mass, center=TRUE, scale=TRUE) +
(1| female) + (1| tank) ,
family=gaussian())
model1.p3 <- brm(f.model1.p3,
data = lat_resp_dat2,
prior = rate.priors2,
warmup = 500,
iter = 5000,
seed=123,
cores=2,
save_pars = save_pars(all=TRUE),
sample_prior = "yes",
chains = 2,
thin = 5,
control = list(adapt_delta=0.95)) ## Compiling Stan program...## Start samplingDHARMa
#--- distribution check ---#
pp_check(model1.p3, type = 'dens_overlay_grouped', ndraws=150, group="region")
#--- DHARMa checks ---#
preds <- posterior_predict(model1.p3, ndraws=250, summary=FALSE)
model1.p3_resids <- createDHARMa(simulatedResponse = t(preds),
observedResponse = lat_resp_dat2$rate.output,
fittedPredictedResponse = apply(preds, 2, mean),
integerResponse = 'student')
plot(model1.p3_resids) ; testDispersion(model1.p3_resids)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1078090, p-value < 2.2e-16
## alternative hypothesis: two.sidedMCMC diagnostics
## 'pars' not specified. Showing first 10 parameters by default.
## 'pars' not specified. Showing first 10 parameters by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## 'pars' not specified. Showing first 10 parameters by default.
LOO
## Output of model 'model1.p1':
##
## Computed from 1800 by 4776 log-likelihood matrix.
##
## Estimate SE
## elpd_loo 1452.7 84.8
## p_loo 60.6 2.7
## looic -2905.4 169.6
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.8, 1.1]).
##
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
##
## Output of model 'model1.p2':
##
## Computed from 1800 by 4776 log-likelihood matrix.
##
## Estimate SE
## elpd_loo 1628.3 88.7
## p_loo 60.9 2.8
## looic -3256.6 177.3
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.6, 1.1]).
##
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
##
## Output of model 'model1.p3':
##
## Computed from 1800 by 4776 log-likelihood matrix.
##
## Estimate SE
## elpd_loo 1701.8 84.3
## p_loo 61.8 2.7
## looic -3403.7 168.7
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.7, 1.1]).
##
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
##
## Model comparisons:
## elpd_diff se_diff
## model1.p3 0.0 0.0
## model1.p2 -73.5 9.2
## model1.p1 -249.1 21.5Fit model [fixed factors - GAM]
GAM model k=4
GAM model k=4
f.model1.gamk4 <- bf(rate.output2 ~ 1 +
t2(dev.temp,region,Temperature, k=4, bs="fs", by=dev.temp.region) +
scale(mass, center=TRUE, scale=TRUE) +
(1| female) + (1| tank),
family=gaussian())
model1.gamk4 <- brm(f.model1.gamk4,
data = lat_resp_dat2,
prior = rate.priors2,
warmup = 500,
iter = 5000,
seed=123,
cores=2,
save_pars = save_pars(all=TRUE),
sample_prior = "yes",
chains = 2,
thin = 5,
control = list(adapt_delta=0.95)) ## Compiling Stan program...## Start sampling## Warning: There were 2 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.## Warning: There were 899 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
## https://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded## Warning: Examine the pairs() plot to diagnose sampling problemsDHARMa
#--- distribution check ---#
pp_check(model1.gamk4, type = 'dens_overlay_grouped', ndraws=150, group="dev.temp.region")
#--- DHARMa checks ---#
preds <- posterior_predict(model1.gamk4, ndraws=250, summary=FALSE)
model1.gamk4_resids <- createDHARMa(simulatedResponse = t(preds),
observedResponse = lat_resp_dat2$rate.output,
fittedPredictedResponse = apply(preds, 2, mean),
integerResponse = 'student')
plot(model1.gamk4_resids) ; testDispersion(model1.gamk4_resids)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1199428, p-value < 2.2e-16
## alternative hypothesis: two.sidedMCMC diagnostics
## 'pars' not specified. Showing first 10 parameters by default.
## 'pars' not specified. Showing first 10 parameters by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## 'pars' not specified. Showing first 10 parameters by default.
### {-}
GAM model k=3
GAM model k=3
f.model1.gamk3 <- bf(rate.output2 ~ 1 +
t2(dev.temp,region,Temperature, k=3, bs="fs", by=dev.temp.region) +
scale(mass, center=TRUE, scale=TRUE) +
(1| female) + (1| tank),
family=gaussian())
model1.gamk3 <- brm(f.model1.gamk3,
data = lat_resp_dat2,
prior = rate.priors2,
warmup = 500,
iter = 5000,
seed=123,
cores=2,
save_pars = save_pars(all=TRUE),
sample_prior = "yes",
chains = 2,
thin = 5,
control = list(adapt_delta=0.95)) ## Compiling Stan program...## Start sampling## Warning: There were 25 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.## Warning: There were 1422 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
## https://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded## Warning: Examine the pairs() plot to diagnose sampling problemsDHARMa
#--- distribution check ---#
pp_check(model1.gamk3, type = 'dens_overlay_grouped', ndraws=150, group="region")
#--- DHARMa checks ---#
preds <- posterior_predict(model1.gamk3, ndraws=250, summary=FALSE)
model1.gamk3_resids <- createDHARMa(simulatedResponse = t(preds),
observedResponse = lat_resp_dat2$rate.output,
fittedPredictedResponse = apply(preds, 2, mean),
integerResponse = 'student')
plot(model1.gamk3_resids) ; testDispersion(model1.gamk3_resids)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1160466, p-value < 2.2e-16
## alternative hypothesis: two.sidedMCMC diagnostics
## 'pars' not specified. Showing first 10 parameters by default.
## 'pars' not specified. Showing first 10 parameters by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## 'pars' not specified. Showing first 10 parameters by default.
### {-}
## Warning in AIC.default(model1.gamk3, model1.gamk4): models are not all fitted
## to the same number of observationsPartial effects plots
Model investigation
summary
## Warning: There were 2 divergent transitions after warmup. Increasing
## adapt_delta above 0.95 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: rate.output2 ~ 1 + t2(dev.temp, region, Temperature, k = 4, bs = "fs", by = dev.temp.region) + scale(mass, center = TRUE, scale = TRUE) + (1 | female) + (1 | tank)
## Data: lat_resp_dat2 (Number of observations: 4776)
## Draws: 2 chains, each with iter = 5000; warmup = 500; thin = 5;
## total post-warmup draws = 1800
##
## Smooth Terms:
## Estimate Est.Error
## sds(t2dev.tempregionTemperaturedev.temp.region28_core_1) 1.49 1.26
## sds(t2dev.tempregionTemperaturedev.temp.region28_core_2) 1.52 1.26
## sds(t2dev.tempregionTemperaturedev.temp.region28_leading_1) 2.33 1.37
## sds(t2dev.tempregionTemperaturedev.temp.region28_leading_2) 2.27 1.33
## sds(t2dev.tempregionTemperaturedev.temp.region30_core_1) 1.94 1.12
## sds(t2dev.tempregionTemperaturedev.temp.region30_core_2) 1.99 1.25
## sds(t2dev.tempregionTemperaturedev.temp.region30_leading_1) 2.17 1.25
## sds(t2dev.tempregionTemperaturedev.temp.region30_leading_2) 2.26 1.27
## sds(t2dev.tempregionTemperaturedev.temp.region31_core_1) 1.42 1.06
## sds(t2dev.tempregionTemperaturedev.temp.region31_core_2) 1.54 1.23
## sds(t2dev.tempregionTemperaturedev.temp.region31_leading_1) 1.70 1.00
## sds(t2dev.tempregionTemperaturedev.temp.region31_leading_2) 1.71 1.21
## l-95% CI u-95% CI
## sds(t2dev.tempregionTemperaturedev.temp.region28_core_1) 0.46 4.10
## sds(t2dev.tempregionTemperaturedev.temp.region28_core_2) 0.43 4.44
## sds(t2dev.tempregionTemperaturedev.temp.region28_leading_1) 0.86 5.70
## sds(t2dev.tempregionTemperaturedev.temp.region28_leading_2) 0.88 5.63
## sds(t2dev.tempregionTemperaturedev.temp.region30_core_1) 0.70 4.77
## sds(t2dev.tempregionTemperaturedev.temp.region30_core_2) 0.75 5.15
## sds(t2dev.tempregionTemperaturedev.temp.region30_leading_1) 0.79 5.35
## sds(t2dev.tempregionTemperaturedev.temp.region30_leading_2) 0.85 5.61
## sds(t2dev.tempregionTemperaturedev.temp.region31_core_1) 0.41 4.36
## sds(t2dev.tempregionTemperaturedev.temp.region31_core_2) 0.46 4.50
## sds(t2dev.tempregionTemperaturedev.temp.region31_leading_1) 0.58 4.19
## sds(t2dev.tempregionTemperaturedev.temp.region31_leading_2) 0.55 4.76
## Rhat Bulk_ESS
## sds(t2dev.tempregionTemperaturedev.temp.region28_core_1) 1.00 1713
## sds(t2dev.tempregionTemperaturedev.temp.region28_core_2) 1.00 1567
## sds(t2dev.tempregionTemperaturedev.temp.region28_leading_1) 1.00 1550
## sds(t2dev.tempregionTemperaturedev.temp.region28_leading_2) 1.00 1819
## sds(t2dev.tempregionTemperaturedev.temp.region30_core_1) 1.00 1616
## sds(t2dev.tempregionTemperaturedev.temp.region30_core_2) 1.00 1855
## sds(t2dev.tempregionTemperaturedev.temp.region30_leading_1) 1.00 1561
## sds(t2dev.tempregionTemperaturedev.temp.region30_leading_2) 1.00 1620
## sds(t2dev.tempregionTemperaturedev.temp.region31_core_1) 1.00 1586
## sds(t2dev.tempregionTemperaturedev.temp.region31_core_2) 1.00 1787
## sds(t2dev.tempregionTemperaturedev.temp.region31_leading_1) 1.00 1663
## sds(t2dev.tempregionTemperaturedev.temp.region31_leading_2) 1.00 1780
## Tail_ESS
## sds(t2dev.tempregionTemperaturedev.temp.region28_core_1) 1631
## sds(t2dev.tempregionTemperaturedev.temp.region28_core_2) 1503
## sds(t2dev.tempregionTemperaturedev.temp.region28_leading_1) 1630
## sds(t2dev.tempregionTemperaturedev.temp.region28_leading_2) 1872
## sds(t2dev.tempregionTemperaturedev.temp.region30_core_1) 1536
## sds(t2dev.tempregionTemperaturedev.temp.region30_core_2) 1709
## sds(t2dev.tempregionTemperaturedev.temp.region30_leading_1) 1610
## sds(t2dev.tempregionTemperaturedev.temp.region30_leading_2) 1687
## sds(t2dev.tempregionTemperaturedev.temp.region31_core_1) 1792
## sds(t2dev.tempregionTemperaturedev.temp.region31_core_2) 1743
## sds(t2dev.tempregionTemperaturedev.temp.region31_leading_1) 1477
## sds(t2dev.tempregionTemperaturedev.temp.region31_leading_2) 1591
##
## Group-Level Effects:
## ~female (Number of levels: 15)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.10 0.04 0.03 0.20 1.00 878 1032
##
## ~tank (Number of levels: 52)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.14 0.02 0.11 0.18 1.00 1499 1454
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat
## Intercept 0.49 0.08 0.32 0.64 1.00
## scalemasscenterEQTRUEscaleEQTRUE 0.26 0.01 0.25 0.27 1.00
## Bulk_ESS Tail_ESS
## Intercept 1635 1644
## scalemasscenterEQTRUEscaleEQTRUE 1750 1717
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.16 0.00 0.16 0.16 1.00 1883 1743
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).gather draws
#model1.re.wo |> get_variables()
model1.gamk4 |> gather_draws(`b_.*|sigma`, regex =TRUE) |>
median_hdci()
emmeans - pariwise
model1.gamk4 |> emmeans(pairwise ~ region*dev.temp, type="response") |> pairs(by="dev.temp") |> summary(infer=TRUE)## NOTE: A nesting structure was detected in the fitted model:
## dev.temp.region %in% (dev.temp*region)model1.gamk4 |> emmeans(pairwise ~ region*dev.temp, type="response") |> pairs(by="region") |> confint()## NOTE: A nesting structure was detected in the fitted model:
## dev.temp.region %in% (dev.temp*region)## NOTE: A nesting structure was detected in the fitted model:
## dev.temp.region %in% (dev.temp*region)## region dev.temp emmean lower.HPD upper.HPD
## core 28 0.531 0.412 0.659
## leading 28 0.550 0.422 0.698
## core 30 0.635 0.526 0.737
## leading 30 0.498 0.352 0.640
## core 31 0.618 0.494 0.725
## leading 31 0.607 0.487 0.730
##
## Point estimate displayed: median
## HPD interval probability: 0.95emtrends - pariwise
## NOTE: A nesting structure was detected in the fitted model:
## dev.temp.region %in% (dev.temp*region)## $emtrends
## dev.temp region Temperature.trend lower.HPD upper.HPD
## 28 core 0.05589 0.0364 0.07286
## 30 core 0.00379 -0.0142 0.01943
## 31 core 0.02995 0.0132 0.04545
## 28 leading -0.01301 -0.0325 0.00544
## 30 leading 0.00435 -0.0153 0.02521
## 31 leading 0.01009 -0.0071 0.02705
##
## Point estimate displayed: median
## HPD interval probability: 0.95
##
## $contrasts
## contrast estimate lower.HPD upper.HPD
## dev.temp28 core - dev.temp30 core 0.052097 0.025961 0.07356
## dev.temp28 core - dev.temp31 core 0.025937 0.000936 0.04901
## dev.temp28 core - dev.temp28 leading 0.068667 0.043783 0.09549
## dev.temp28 core - dev.temp30 leading 0.051514 0.025472 0.07994
## dev.temp28 core - dev.temp31 leading 0.046443 0.020681 0.06957
## dev.temp30 core - dev.temp31 core -0.026011 -0.049167 -0.00357
## dev.temp30 core - dev.temp28 leading 0.017010 -0.007222 0.04247
## dev.temp30 core - dev.temp30 leading -0.000309 -0.025766 0.02605
## dev.temp30 core - dev.temp31 leading -0.005979 -0.031522 0.01642
## dev.temp31 core - dev.temp28 leading 0.042749 0.019082 0.06782
## dev.temp31 core - dev.temp30 leading 0.025722 -0.000808 0.05269
## dev.temp31 core - dev.temp31 leading 0.019962 -0.004044 0.04130
## dev.temp28 leading - dev.temp30 leading -0.017101 -0.046804 0.00952
## dev.temp28 leading - dev.temp31 leading -0.023032 -0.049326 0.00121
## dev.temp30 leading - dev.temp31 leading -0.005859 -0.032276 0.02083
##
## Point estimate displayed: median
## HPD interval probability: 0.95probabilities
## NOTE: A nesting structure was detected in the fitted model:
## dev.temp.region %in% (dev.temp*region)mtsqrt2 <- mtsqst$contrasts |> gather_emmeans_draws()
mtsqrt2 %>% group_by(contrast) %>% dplyr::summarise(Prob = sum(.value>0)/n())values
values.df <- lat_resp_dat2 |>
group_by(region, dev.temp,dev.temp.region) |>
data_grid(Temperature=seq(from=min(lat_resp_dat2$Temperature),
to=max(lat_resp_dat2$Temperature),
length.out=100),
mass=0.00148299) |>
add_fitted_draws(model1.gamk4,
n=300,
re_formula=NA) |>
filter(Temperature < 34) |>
group_by(dev.temp.region) |>
summarise(oxygen_uptake_values = median(.value)) |>
ungroup() |>
separate(dev.temp.region, c("temperature","region"), remove=FALSE); values.df## Warning: `fitted_draws` and `add_fitted_draws` are deprecated as their names were confusing.
## - Use [add_]epred_draws() to get the expectation of the posterior predictive.
## - Use [add_]linpred_draws() to get the distribution of the linear predictor.
## - For example, you used [add_]fitted_draws(..., scale = "response"), which
## means you most likely want [add_]epred_draws(...).
## NOTE: When updating to the new functions, note that the `model` parameter is now
## named `object` and the `n` parameter is now named `ndraws`.summary figure
resp.plot <- lat_resp_dat2 |>
group_by(region, dev.temp,dev.temp.region) |>
data_grid(Temperature=seq(from=min(lat_resp_dat2$Temperature),
to=max(lat_resp_dat2$Temperature),
length.out=100),
mass=0.00148299) |>
add_fitted_draws(model1.gamk4,
n=300,
re_formula=NA) |>
unite("exp_group", c("region","dev.temp"), remove=FALSE) |>
ggplot(aes(x=Temperature, color=exp_group)) +
geom_line(aes(y=.value, group=paste(exp_group, .draw)), alpha=0.05) +
geom_smooth(aes(y=.value, color=exp_group), se=FALSE, size=1.5) +
theme_classic() +
facet_wrap(~dev.temp); resp.plot
resp.plot2 <- lat_resp_dat2 |>
group_by(region, dev.temp,dev.temp.region) |>
data_grid(Temperature=seq(from=min(lat_resp_dat2$Temperature),
to=max(lat_resp_dat2$Temperature),
length.out=100),
mass=0.00148299) |>
add_fitted_draws(model1.gamk4,
n=300,
re_formula=NA) |>
unite("exp_group", c("region","dev.temp"), remove=FALSE) |>
ggplot(aes(x=Temperature, color=exp_group)) +
geom_line(aes(y=.value, group=paste(exp_group, .draw)), alpha=0.05) +
geom_smooth(aes(y=.value, color=exp_group), se=FALSE, size=1.5) +
theme_classic() +
facet_wrap(~region); resp.plot2
resp.plot3 <- lat_resp_dat2 |>
group_by(region, dev.temp,dev.temp.region) |>
data_grid(Temperature=seq(from=min(lat_resp_dat2$Temperature),
to=max(lat_resp_dat2$Temperature),
length.out=100),
mass=0.00148299) |>
add_fitted_draws(model1.gamk4,
n=300,
re_formula=NA) |>
unite("exp_group", c("region","dev.temp"), remove=FALSE) |>
ggplot(aes(x=Temperature, color=exp_group)) +
geom_line(aes(y=.value, group=paste(exp_group, .draw)), alpha=0.05) +
geom_smooth(aes(y=.value, color=exp_group), se=FALSE, size=1.5) +
theme_classic(); resp.plot3
Pre-inflection point
create data frame
Split posterior of draws from main model into data prior to 33.5 C (earliest inflection point)
priors
pre_infl |>
group_by(dev.temp.region) |>
summarise(mean(na.omit(rate.output2)),
sd(na.omit(rate.output2)))linear model
f.pre.model1 <- bf(rate.output2 ~ 1 +
dev.temp*region*Temperature +
scale(mass, center=TRUE, scale=TRUE) +
(1| female) + (1| tank),
family=gaussian())
pre.model1 <- brm(f.pre.model1,
data = pre_infl,
prior = pre.priors,
warmup = 500,
iter = 5000,
seed=123,
cores=2,
save_pars = save_pars(all=TRUE),
sample_prior = "yes",
chains = 2,
thin = 5,
control = list(adapt_delta=0.95)) ## Compiling Stan program...## Start sampling## Warning: There were 758 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
## https://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded## Warning: Examine the pairs() plot to diagnose sampling problemsModel investigation
summary
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: rate.output2 ~ 1 + dev.temp * region * Temperature + scale(mass, center = TRUE, scale = TRUE) + (1 | female) + (1 | tank)
## Data: pre_infl (Number of observations: 2439)
## Draws: 2 chains, each with iter = 5000; warmup = 500; thin = 5;
## total post-warmup draws = 1800
##
## Group-Level Effects:
## ~female (Number of levels: 15)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.10 0.04 0.03 0.20 1.00 1413 1566
##
## ~tank (Number of levels: 52)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.14 0.02 0.11 0.19 1.00 1410 1336
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat
## Intercept 0.30 0.14 0.03 0.60 1.00
## dev.temp30 -0.18 0.17 -0.51 0.14 1.00
## dev.temp31 0.13 0.16 -0.17 0.46 1.00
## regionleading 0.12 0.16 -0.21 0.42 1.00
## Temperature 0.01 0.00 -0.00 0.01 1.00
## scalemasscenterEQTRUEscaleEQTRUE 0.25 0.01 0.24 0.26 1.00
## dev.temp30:regionleading 0.14 0.20 -0.25 0.53 1.00
## dev.temp31:regionleading 0.14 0.19 -0.23 0.49 1.00
## dev.temp30:Temperature 0.01 0.01 0.00 0.02 1.00
## dev.temp31:Temperature -0.00 0.01 -0.01 0.01 1.00
## regionleading:Temperature -0.00 0.01 -0.01 0.01 1.00
## dev.temp30:regionleading:Temperature -0.01 0.01 -0.03 0.00 1.00
## dev.temp31:regionleading:Temperature -0.01 0.01 -0.02 0.01 1.00
## Bulk_ESS Tail_ESS
## Intercept 1993 1737
## dev.temp30 1959 1738
## dev.temp31 1846 1679
## regionleading 1615 1795
## Temperature 1846 1743
## scalemasscenterEQTRUEscaleEQTRUE 1794 1667
## dev.temp30:regionleading 1698 1708
## dev.temp31:regionleading 1776 1635
## dev.temp30:Temperature 1941 1747
## dev.temp31:Temperature 1769 1664
## regionleading:Temperature 1565 1648
## dev.temp30:regionleading:Temperature 1723 1683
## dev.temp31:regionleading:Temperature 1738 1678
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.12 0.00 0.12 0.12 1.00 1666 1701
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).gather draws
#model1.re.wo |> get_variables()
pre.model1 |> gather_draws(`b_.*|sigma`, regex =TRUE) |>
median_hdci()
emmeans - pariwise
pre.model1 |> emmeans(pairwise ~ region*dev.temp, type="response") |> pairs(by="dev.temp") |> summary(infer=TRUE)## NOTE: Results may be misleading due to involvement in interactionspre.model1 |> emmeans(pairwise ~ region*dev.temp, type="response") |> pairs(by="region") |> confint()## NOTE: Results may be misleading due to involvement in interactions## NOTE: Results may be misleading due to involvement in interactions## region dev.temp emmean lower.HPD upper.HPD
## core 28 0.500 0.370 0.617
## leading 28 0.563 0.429 0.705
## core 30 0.650 0.542 0.760
## leading 30 0.471 0.343 0.626
## core 31 0.597 0.473 0.708
## leading 31 0.599 0.462 0.724
##
## Point estimate displayed: median
## HPD interval probability: 0.95Summary plots
resp.plot <- lat_resp_dat2 |>
group_by(region, dev.temp,dev.temp.region) |>
data_grid(Temperature=seq(from=min(lat_resp_dat2$Temperature),
to=33.5,
length.out=100),
mass=0.00148299) |>
add_fitted_draws(pre.model1,
n=500,
re_formula=NA) |>
unite("exp_group", c("region","dev.temp"), remove=FALSE) |>
ggplot(aes(x=Temperature, color=exp_group)) +
geom_line(aes(y=.value, group=paste(exp_group, .draw)), alpha=0.05) +
geom_smooth(aes(y=.value, color=exp_group), se=FALSE, size=1.5) +
theme_classic(); resp.plot
Post-inflection point analysis
priors
post_infl |>
group_by(dev.temp.region) |>
summarise(mean(na.omit(rate.output2)),
sd(na.omit(rate.output2)))model
f.post.model1 <- bf(rate.output2 ~ 1 +
dev.temp*region*Temperature +
scale(mass, center=TRUE, scale=TRUE) +
(1| female) + (1| tank),
family=gaussian())
post.model1 <- brm(f.post.model1,
data = post_infl,
prior = post.priors,
warmup = 500,
iter = 5000,
seed=123,
cores=2,
save_pars = save_pars(all=TRUE),
sample_prior = "yes",
chains = 2,
thin = 5,
control = list(adapt_delta=0.95)) ## Compiling Stan program...## Start samplingModel investigation
summary
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: rate.output2 ~ 1 + dev.temp * region * Temperature + scale(mass, center = TRUE, scale = TRUE) + (1 | female) + (1 | tank)
## Data: post_infl (Number of observations: 1475)
## Draws: 2 chains, each with iter = 5000; warmup = 500; thin = 5;
## total post-warmup draws = 1800
##
## Group-Level Effects:
## ~female (Number of levels: 15)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.09 0.05 0.01 0.20 1.00 1111 1422
##
## ~tank (Number of levels: 52)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.17 0.02 0.13 0.22 1.00 1594 1721
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat
## Intercept -7.17 0.32 -7.80 -6.56 1.00
## dev.temp30 -0.17 0.19 -0.54 0.21 1.00
## dev.temp31 0.23 0.19 -0.14 0.59 1.00
## regionleading -0.18 0.19 -0.56 0.20 1.00
## Temperature 0.22 0.01 0.21 0.24 1.00
## scalemasscenterEQTRUEscaleEQTRUE 0.28 0.01 0.26 0.30 1.00
## dev.temp30:regionleading -0.09 0.19 -0.45 0.28 1.00
## dev.temp31:regionleading 0.02 0.20 -0.36 0.43 1.00
## dev.temp30:Temperature 0.01 0.01 -0.01 0.02 1.00
## dev.temp31:Temperature -0.01 0.01 -0.02 0.01 1.00
## regionleading:Temperature 0.01 0.01 -0.01 0.02 1.00
## dev.temp30:regionleading:Temperature -0.00 0.01 -0.01 0.01 1.00
## dev.temp31:regionleading:Temperature -0.00 0.01 -0.01 0.01 1.00
## Bulk_ESS Tail_ESS
## Intercept 1638 1785
## dev.temp30 1710 1646
## dev.temp31 1883 1708
## regionleading 1735 1697
## Temperature 1627 1610
## scalemasscenterEQTRUEscaleEQTRUE 1767 1685
## dev.temp30:regionleading 1624 1708
## dev.temp31:regionleading 1907 1808
## dev.temp30:Temperature 1688 1598
## dev.temp31:Temperature 1901 1613
## regionleading:Temperature 1699 1716
## dev.temp30:regionleading:Temperature 1443 1526
## dev.temp31:regionleading:Temperature 1876 1704
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.20 0.00 0.20 0.21 1.00 1592 1405
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).gather draws
#model1.re.wo |> get_variables()
post.model1 |> gather_draws(`b_.*|sigma`, regex =TRUE) |>
median_hdci()
emmeans - pariwise
## $emtrends
## dev.temp region Temperature.trend lower.HPD upper.HPD
## 28 core 0.222 0.205 0.240
## 30 core 0.228 0.210 0.247
## 31 core 0.217 0.200 0.236
## 28 leading 0.228 0.211 0.246
## 30 leading 0.235 0.212 0.254
## 31 leading 0.222 0.202 0.241
##
## Point estimate displayed: median
## HPD interval probability: 0.95
##
## $contrasts
## contrast estimate lower.HPD upper.HPD
## dev.temp28 core - dev.temp30 core -0.006634 -0.01665 0.00580
## dev.temp28 core - dev.temp31 core 0.005451 -0.00620 0.01579
## dev.temp28 core - dev.temp28 leading -0.006059 -0.01806 0.00556
## dev.temp28 core - dev.temp30 leading -0.012625 -0.03053 0.00557
## dev.temp28 core - dev.temp31 leading 0.000468 -0.01708 0.01853
## dev.temp30 core - dev.temp31 core 0.012027 -0.00216 0.02650
## dev.temp30 core - dev.temp28 leading 0.000288 -0.01667 0.01554
## dev.temp30 core - dev.temp30 leading -0.006018 -0.02140 0.00908
## dev.temp30 core - dev.temp31 leading 0.006744 -0.01352 0.02615
## dev.temp31 core - dev.temp28 leading -0.011820 -0.02571 0.00631
## dev.temp31 core - dev.temp30 leading -0.017863 -0.04033 0.00123
## dev.temp31 core - dev.temp31 leading -0.005220 -0.02090 0.00978
## dev.temp28 leading - dev.temp30 leading -0.006391 -0.02154 0.00864
## dev.temp28 leading - dev.temp31 leading 0.006235 -0.00854 0.02151
## dev.temp30 leading - dev.temp31 leading 0.012649 -0.00721 0.03365
##
## Point estimate displayed: median
## HPD interval probability: 0.95Summary plots
resp.plot <- lat_resp_dat2 |>
group_by(region, dev.temp,dev.temp.region) |>
data_grid(Temperature=seq(from=34.6,
to=max(lat_resp_dat2$Temperature),
length.out=100),
mass=0.00148299) |>
add_fitted_draws(post.model1,
n=500,
re_formula=NA) |>
unite("exp_group", c("region","dev.temp"), remove=FALSE) |>
ggplot(aes(x=Temperature, color=exp_group)) +
geom_line(aes(y=.value, group=paste(exp_group, .draw)), alpha=0.05) +
geom_smooth(aes(y=.value, color=exp_group), se=FALSE, size=1.5) +
theme_classic(); resp.plot
Piecewise linear regression
Low-latitude - 28.5C
posterior_draws <- lat_resp_dat2 |>
group_by(region, dev.temp,dev.temp.region) |>
data_grid(Temperature=seq(from=min(lat_resp_dat2$Temperature),
to=max(lat_resp_dat2$Temperature),
length.out=100),
mass=0.00148299) |>
add_linpred_draws(model1.gamk4,
n=300,
re_formula=NA) |>
unite("exp_group", c("region","dev.temp"), remove=FALSE) |>
filter(dev.temp.region == "28_core")
spline_fit <- smooth.spline(posterior_draws$Temperature, posterior_draws$.linpred, nknots = 4)
ggplot(posterior_draws, aes(x=Temperature, y=.linpred)) +
geom_point(color="black", alpha=0.05, size=1) +
geom_smooth(data=as.data.frame(spline_fit$data), aes(x=x, y=y,
color="red"),
size=2) ## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
strucchange analysis
spline.data <- as.data.frame(x= spline_fit$x) |>
mutate(y = spline_fit$y)|>
rename(x = `spline_fit$x`)
plot(spline.data)
##
## Optimal (m+1)-segment partition:
##
## Call:
## breakpoints.formula(formula = spline.data$y ~ 1)
##
## Breakpoints at observation number:
##
## m = 1 73
## m = 2 64 83
## m = 3 55 70 85
## m = 4 40 55 70 85
## m = 5 15 40 55 70 85
##
## Corresponding to breakdates:
##
## m = 1 0.73
## m = 2 0.64 0.83
## m = 3 0.55 0.7 0.85
## m = 4 0.4 0.55 0.7 0.85
## m = 5 0.15 0.4 0.55 0.7 0.85
##
## Fit:
##
## m 0 1 2 3 4 5
## RSS 2.7785 0.5048 0.2062 0.1216 0.1174 0.1116
## BIC -65.3272 -226.6618 -306.9730 -350.5902 -344.8844 -340.7792
## [1] 33.86963 35.04506 36.14703segement analysis
Use out from strucchange summary to number of breakpoints as well as approximate location of breakpoints
## [1] 33.86963 35.04506 36.14703my.seg <- segmented(my.lm,
seg.Z = ~x,
psi=list(x=spline.data[breakpoint.values, c("x")]))
#summary(my.seg)
my.seg$psi ; slope(my.seg)## Initial Est. St.Err
## psi1.x 33.86963 32.35335 0.03176730
## psi2.x 35.04506 33.76564 0.02974222
## psi3.x 36.14703 35.09270 0.02693903## $x
## Est. St.Err. t value CI(95%).l CI(95%).u
## slope1 -0.019400 0.00076011 -25.522 -0.020909 -0.017890
## slope2 0.036286 0.00174150 20.836 0.032828 0.039745
## slope3 0.103900 0.00188890 55.003 0.100140 0.107650
## slope4 0.171910 0.00087702 196.010 0.170160 0.173650my.fitted <- fitted(my.seg)
my.model <- data.frame(temperature = spline.data$x, rate.output= my.fitted)
ggplot(my.model, aes(x=temperature, y=rate.output)) + geom_line()
ggplot(posterior_draws, aes(x=Temperature, y=.linpred, color=.linpred)) +
geom_point(alpha=0.05, size=1) +
geom_smooth(data=as.data.frame(spline_fit$data), aes(x=x, y=y), color="red",
size=3) + ylab("rate.output")+
geom_vline(xintercept = my.seg$psi[1,2]) +
geom_vline(xintercept = my.seg$psi[2,2])+
geom_vline(xintercept = my.seg$psi[3,2])+
geom_vline(xintercept = my.seg$psi[1,1], color="blue") +
geom_vline(xintercept = my.seg$psi[2,1], color="blue") +
geom_vline(xintercept = my.seg$psi[3,1], color="blue") +
geom_line(data=my.model, aes(x=temperature, y=rate.output), color="black", size=2) +
theme_classic()## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
Low-latitude - 30C
posterior_draws <- lat_resp_dat2 |>
group_by(region, dev.temp,dev.temp.region) |>
data_grid(Temperature=seq(from=min(lat_resp_dat2$Temperature),
to=max(lat_resp_dat2$Temperature),
length.out=100),
mass=0.00148299) |>
add_linpred_draws(model1.gamk4,
n=300,
re_formula=NA) |>
unite("exp_group", c("region","dev.temp"), remove=FALSE) |>
filter(dev.temp.region == "30_core")
spline_fit <- smooth.spline(posterior_draws$Temperature, posterior_draws$.linpred, nknots = 4)
ggplot(posterior_draws, aes(x=Temperature, y=.linpred)) +
geom_point(color="black", alpha=0.05, size=1) +
geom_smooth(data=as.data.frame(spline_fit$data), aes(x=x, y=y,
color="red"),
size=2) ## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
strucchange analysis
spline.data <- as.data.frame(x= spline_fit$x) |>
mutate(y = spline_fit$y)|>
rename(x = `spline_fit$x`)
plot(spline.data)
##
## Optimal (m+1)-segment partition:
##
## Call:
## breakpoints.formula(formula = spline.data$y ~ 1)
##
## Breakpoints at observation number:
##
## m = 1 78
## m = 2 70 85
## m = 3 18 70 85
## m = 4 15 55 70 85
## m = 5 15 30 55 70 85
##
## Corresponding to breakdates:
##
## m = 1 0.78
## m = 2 0.7 0.85
## m = 3 0.18 0.7 0.85
## m = 4 0.15 0.55 0.7 0.85
## m = 5 0.15 0.3 0.55 0.7 0.85
##
## Fit:
##
## m 0 1 2 3 4 5
## RSS 2.9454 0.5269 0.2322 0.1865 0.1643 0.1634
## BIC -59.4938 -222.3828 -295.1253 -307.8064 -311.2763 -302.6371
## [1] 33.86963 35.04506 36.14703segement analysis
Use out from strucchange summary to number of breakpoints as well as approximate location of breakpoints
## [1] 33.86963 35.04506 36.14703my.seg <- segmented(my.lm,
seg.Z = ~x,
psi=list(x=spline.data[breakpoint.values, c("x")]))
#summary(my.seg)
my.seg$psi ; slope(my.seg)## Initial Est. St.Err
## psi1.x 33.86963 34.42398 0.07889169
## psi2.x 35.04506 35.07879 0.12326095
## psi3.x 36.14703 35.83758 0.14302773## $x
## Est. St.Err. t value CI(95%).l CI(95%).u
## slope1 0.017059 0.00099279 17.1830 0.015088 0.019031
## slope2 0.107630 0.01762500 6.1065 0.072621 0.142630
## slope3 0.182670 0.01503100 12.1530 0.152810 0.212520
## slope4 0.238360 0.00529410 45.0230 0.227840 0.248870my.fitted <- fitted(my.seg)
my.model <- data.frame(temperature = spline.data$x, rate.output= my.fitted)
ggplot(my.model, aes(x=temperature, y=rate.output)) + geom_line()
ggplot(posterior_draws, aes(x=Temperature, y=.linpred, color=.linpred)) +
geom_point(alpha=0.05, size=1) +
geom_smooth(data=as.data.frame(spline_fit$data), aes(x=x, y=y), color="red",
size=3) + ylab("rate.output")+
geom_vline(xintercept = my.seg$psi[1,2]) +
geom_vline(xintercept = my.seg$psi[2,2])+
geom_vline(xintercept = my.seg$psi[3,2])+
geom_vline(xintercept = my.seg$psi[1,1], color="blue") +
geom_vline(xintercept = my.seg$psi[2,1], color="blue") +
geom_vline(xintercept = my.seg$psi[3,1], color="blue") +
geom_line(data=my.model, aes(x=temperature, y=rate.output), color="black", size=2) +
theme_classic() ## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
Low-latitude - 31C
posterior_draws <- lat_resp_dat2 |>
group_by(region, dev.temp,dev.temp.region) |>
data_grid(Temperature=seq(from=min(lat_resp_dat2$Temperature),
to=max(lat_resp_dat2$Temperature),
length.out=100),
mass=0.00148299) |>
add_linpred_draws(model1.gamk4,
n=300,
re_formula=NA) |>
unite("exp_group", c("region","dev.temp"), remove=FALSE) |>
filter(dev.temp.region == "31_core")
spline_fit <- smooth.spline(posterior_draws$Temperature, posterior_draws$.linpred, nknots = 4)
ggplot(posterior_draws, aes(x=Temperature, y=.linpred)) +
geom_point(color="black", alpha=0.05, size=1) +
geom_smooth(data=as.data.frame(spline_fit$data), aes(x=x, y=y,
color="red"),
size=2) ## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
strucchange analysis
spline.data <- as.data.frame(x= spline_fit$x) |>
mutate(y = spline_fit$y)|>
rename(x = `spline_fit$x`)
plot(spline.data)
##
## Optimal (m+1)-segment partition:
##
## Call:
## breakpoints.formula(formula = spline.data$y ~ 1)
##
## Breakpoints at observation number:
##
## m = 1 75
## m = 2 66 84
## m = 3 55 70 85
## m = 4 40 55 70 85
## m = 5 15 40 55 70 85
##
## Corresponding to breakdates:
##
## m = 1 0.75
## m = 2 0.66 0.84
## m = 3 0.55 0.7 0.85
## m = 4 0.4 0.55 0.7 0.85
## m = 5 0.15 0.4 0.55 0.7 0.85
##
## Fit:
##
## m 0 1 2 3 4 5
## RSS 1.90040 0.34415 0.14201 0.08957 0.08432 0.08360
## BIC -103.31265 -264.97683 -344.28312 -381.16143 -377.99555 -369.63983
## [1] 33.94309 35.04506 36.14703segement analysis
Use out from strucchange summary to number of breakpoints as well as approximate location of breakpoints
## [1] 33.94309 35.04506 36.14703my.seg <- segmented(my.lm,
seg.Z = ~x,
psi=list(x=spline.data[breakpoint.values, c("x")]))
#summary(my.seg)
my.seg$psi ; slope(my.seg)## Initial Est. St.Err
## psi1.x 33.94309 33.32652 0.02649379
## psi2.x 35.04506 34.50521 0.02054399
## psi3.x 36.14703 35.56366 0.02483205## $x
## Est. St.Err. t value CI(95%).l CI(95%).u
## slope1 0.0065759 0.00028252 23.276 0.0060148 0.007137
## slope2 0.0473010 0.00137960 34.286 0.0445610 0.050041
## slope3 0.1127400 0.00152030 74.157 0.1097200 0.115760
## slope4 0.1620300 0.00079966 202.630 0.1604400 0.163620my.fitted <- fitted(my.seg)
my.model <- data.frame(temperature = spline.data$x, rate.output= my.fitted)
ggplot(my.model, aes(x=temperature, y=rate.output)) + geom_line()
ggplot(posterior_draws, aes(x=Temperature, y=.linpred, color=.linpred)) +
geom_point(alpha=0.05, size=1) +
geom_smooth(data=as.data.frame(spline_fit$data), aes(x=x, y=y), color="red",
size=3) + ylab("rate.output")+
geom_vline(xintercept = my.seg$psi[1,2]) +
geom_vline(xintercept = my.seg$psi[2,2])+
geom_vline(xintercept = my.seg$psi[3,2])+
geom_vline(xintercept = my.seg$psi[1,1], color="blue") +
geom_vline(xintercept = my.seg$psi[2,1], color="blue") +
geom_vline(xintercept = my.seg$psi[3,1], color="blue") +
geom_line(data=my.model, aes(x=temperature, y=rate.output), color="black", size=2) +
theme_classic()## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
High-latitude - 28C
posterior_draws <- lat_resp_dat2 |>
group_by(region, dev.temp,dev.temp.region) |>
data_grid(Temperature=seq(from=min(lat_resp_dat2$Temperature),
to=max(lat_resp_dat2$Temperature),
length.out=100),
mass=0.00148299) |>
add_linpred_draws(model1.gamk4,
n=300,
re_formula=NA) |>
unite("exp_group", c("region","dev.temp"), remove=FALSE) |>
filter(dev.temp.region == "28_leading")
spline_fit <- smooth.spline(posterior_draws$Temperature, posterior_draws$.linpred, nknots = 4)
ggplot(posterior_draws, aes(x=Temperature, y=.linpred)) +
geom_point(color="black", alpha=0.05, size=1) +
geom_smooth(data=as.data.frame(spline_fit$data), aes(x=x, y=y,
color="red"),
size=2) ## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
strucchange analysis
spline.data <- as.data.frame(x= spline_fit$x) |>
mutate(y = spline_fit$y)|>
rename(x = `spline_fit$x`)
plot(spline.data)
##
## Optimal (m+1)-segment partition:
##
## Call:
## breakpoints.formula(formula = spline.data$y ~ 1)
##
## Breakpoints at observation number:
##
## m = 1 79
## m = 2 70 85
## m = 3 15 70 85
## m = 4 15 55 70 85
## m = 5 15 40 55 70 85
##
## Corresponding to breakdates:
##
## m = 1 0.79
## m = 2 0.7 0.85
## m = 3 0.15 0.7 0.85
## m = 4 0.15 0.55 0.7 0.85
## m = 5 0.15 0.4 0.55 0.7 0.85
##
## Fit:
##
## m 0 1 2 3 4 5
## RSS 3.5327 0.5780 0.2508 0.2388 0.2338 0.2296
## BIC -41.3137 -213.1309 -287.3916 -283.1020 -276.0267 -268.5972
## [1] 33.94309 35.04506 36.14703segement analysis
Use out from strucchange summary to number of breakpoints as well as approximate location of breakpoints
## [1] 33.94309 35.04506 36.14703my.seg <- segmented(my.lm,
seg.Z = ~x,
psi=list(x=spline.data[breakpoint.values, c("x")]))
#summary(my.seg)
my.seg$psi ; slope(my.seg)## Initial Est. St.Err
## psi1.x 33.94309 34.42427 0.07892526
## psi2.x 35.04506 35.07879 0.12333188
## psi3.x 36.14703 35.83765 0.14306228## $x
## Est. St.Err. t value CI(95%).l CI(95%).u
## slope1 0.0034725 0.0012763 2.7208 0.00093769 0.0060072
## slope2 0.1197800 0.0226570 5.2867 0.07478400 0.1647800
## slope3 0.2162000 0.0193220 11.1890 0.17782000 0.2545700
## slope4 0.2877800 0.0068057 42.2850 0.27426000 0.3012900my.fitted <- fitted(my.seg)
my.model <- data.frame(temperature = spline.data$x, rate.output= my.fitted)
ggplot(my.model, aes(x=temperature, y=rate.output)) + geom_line()
ggplot(posterior_draws, aes(x=Temperature, y=.linpred, color=.linpred)) +
geom_point(alpha=0.05, size=1) +
geom_smooth(data=as.data.frame(spline_fit$data), aes(x=x, y=y), color="red",
size=3) + ylab("rate.output")+
geom_vline(xintercept = my.seg$psi[1,2]) +
geom_vline(xintercept = my.seg$psi[2,2])+
geom_vline(xintercept = my.seg$psi[3,2])+
geom_vline(xintercept = my.seg$psi[1,1], color="blue") +
geom_vline(xintercept = my.seg$psi[2,1], color="blue") +
geom_vline(xintercept = my.seg$psi[3,1], color="blue") +
geom_line(data=my.model, aes(x=temperature, y=rate.output), color="black", size=2) +
theme_classic()## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
High-latitude - 30C
posterior_draws <- lat_resp_dat2 |>
group_by(region, dev.temp,dev.temp.region) |>
data_grid(Temperature=seq(from=min(lat_resp_dat2$Temperature),
to=max(lat_resp_dat2$Temperature),
length.out=100),
mass=0.00148299) |>
add_linpred_draws(model1.gamk4,
n=300,
re_formula=NA) |>
unite("exp_group", c("region","dev.temp"), remove=FALSE) |>
filter(dev.temp.region == "30_leading")
spline_fit <- smooth.spline(posterior_draws$Temperature, posterior_draws$.linpred, nknots = 4)
ggplot(posterior_draws, aes(x=Temperature, y=.linpred)) +
geom_point(color="black", alpha=0.05, size=1) +
geom_smooth(data=as.data.frame(spline_fit$data), aes(x=x, y=y,
color="red"),
size=2) ## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
strucchange analysis
spline.data <- as.data.frame(x= spline_fit$x) |>
mutate(y = spline_fit$y)|>
rename(x = `spline_fit$x`)
plot(spline.data)
##
## Optimal (m+1)-segment partition:
##
## Call:
## breakpoints.formula(formula = spline.data$y ~ 1)
##
## Breakpoints at observation number:
##
## m = 1 78
## m = 2 70 85
## m = 3 55 70 85
## m = 4 15 55 70 85
## m = 5 15 37 55 70 85
##
## Corresponding to breakdates:
##
## m = 1 0.78
## m = 2 0.7 0.85
## m = 3 0.55 0.7 0.85
## m = 4 0.15 0.55 0.7 0.85
## m = 5 0.15 0.37 0.55 0.7 0.85
##
## Fit:
##
## m 0 1 2 3 4 5
## RSS 3.8094 0.6279 0.2542 0.2257 0.2236 0.2225
## BIC -33.7721 -204.8428 -286.0614 -288.7430 -280.4479 -271.7575
## [1] 33.94309 35.04506 36.14703segement analysis
Use out from strucchange summary to number of breakpoints as well as approximate location of breakpoints
## [1] 33.94309 35.04506 36.14703my.seg <- segmented(my.lm,
seg.Z = ~x,
psi=list(x=spline.data[breakpoint.values, c("x")]))
#summary(my.seg)
my.seg$psi ; slope(my.seg)## Initial Est. St.Err
## psi1.x 33.94309 34.13055 0.04832474
## psi2.x 35.04506 34.92601 0.05935603
## psi3.x 36.14703 35.75668 0.07378457## $x
## Est. St.Err. t value CI(95%).l CI(95%).u
## slope1 0.002638 0.00079262 3.3282 0.0010638 0.0042122
## slope2 0.100910 0.00938690 10.7500 0.0822690 0.1195600
## slope3 0.205170 0.00938690 21.8570 0.1865300 0.2238100
## slope4 0.279220 0.00354790 78.6990 0.2721700 0.2862700my.fitted <- fitted(my.seg)
my.model <- data.frame(temperature = spline.data$x, rate.output= my.fitted)
ggplot(my.model, aes(x=temperature, y=rate.output)) + geom_line()
ggplot(posterior_draws, aes(x=Temperature, y=.linpred, color=.linpred)) +
geom_point(alpha=0.05, size=1) +
geom_smooth(data=as.data.frame(spline_fit$data), aes(x=x, y=y), color="red",
size=3) + ylab("rate.output")+
geom_vline(xintercept = my.seg$psi[1,2]) +
geom_vline(xintercept = my.seg$psi[2,2])+
geom_vline(xintercept = my.seg$psi[3,2])+
geom_vline(xintercept = my.seg$psi[1,1], color="blue") +
geom_vline(xintercept = my.seg$psi[2,1], color="blue") +
geom_vline(xintercept = my.seg$psi[3,1], color="blue") +
geom_line(data=my.model, aes(x=temperature, y=rate.output), color="black", size=2) +
theme_classic()## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
High-latitude - 31C
posterior_draws <- lat_resp_dat2 |>
group_by(region, dev.temp,dev.temp.region) |>
data_grid(Temperature=seq(from=min(lat_resp_dat2$Temperature),
to=max(lat_resp_dat2$Temperature),
length.out=100),
mass=0.00148299) |>
add_linpred_draws(model1.gamk4,
n=300,
re_formula=NA) |>
unite("exp_group", c("region","dev.temp"), remove=FALSE) |>
filter(dev.temp.region == "31_leading")
spline_fit <- smooth.spline(posterior_draws$Temperature, posterior_draws$.linpred, nknots = 4)
ggplot(posterior_draws, aes(x=Temperature, y=.linpred)) +
geom_point(color="black", alpha=0.05, size=1) +
geom_smooth(data=as.data.frame(spline_fit$data), aes(x=x, y=y,
color="red"),
size=2) ## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
strucchange analysis
spline.data <- as.data.frame(x= spline_fit$x) |>
mutate(y = spline_fit$y)|>
rename(x = `spline_fit$x`)
plot(spline.data)
##
## Optimal (m+1)-segment partition:
##
## Call:
## breakpoints.formula(formula = spline.data$y ~ 1)
##
## Breakpoints at observation number:
##
## m = 1 78
## m = 2 70 85
## m = 3 55 70 85
## m = 4 32 55 70 85
## m = 5 15 32 55 70 85
##
## Corresponding to breakdates:
##
## m = 1 0.78
## m = 2 0.7 0.85
## m = 3 0.55 0.7 0.85
## m = 4 0.32 0.55 0.7 0.85
## m = 5 0.15 0.32 0.55 0.7 0.85
##
## Fit:
##
## m 0 1 2 3 4 5
## RSS 2.1641 0.3569 0.1421 0.1233 0.1224 0.1223
## BIC -90.3197 -261.3436 -344.1864 -349.1651 -340.7615 -331.5675
## [1] 33.94309 35.04506 36.14703segement analysis
Use out from strucchange summary to number of breakpoints as well as approximate location of breakpoints
## [1] 33.94309 35.04506 36.14703my.seg <- segmented(my.lm,
seg.Z = ~x,
psi=list(x=spline.data[breakpoint.values, c("x")]))
#summary(my.seg)
my.seg$psi ; slope(my.seg)## Initial Est. St.Err
## psi1.x 33.94309 33.83815 0.03068568
## psi2.x 35.04506 34.78840 0.02982365
## psi3.x 36.14703 35.72150 0.03962484## $x
## Est. St.Err. t value CI(95%).l CI(95%).u
## slope1 -0.0032692 0.00039822 -8.2096 -0.0040601 -0.0024783
## slope2 0.0621610 0.00328720 18.9100 0.0556320 0.0686900
## slope3 0.1480500 0.00328720 45.0370 0.1415200 0.1545800
## slope4 0.2054600 0.00159820 128.5600 0.2022800 0.2086300my.fitted <- fitted(my.seg)
my.model <- data.frame(temperature = spline.data$x, rate.output= my.fitted)
ggplot(my.model, aes(x=temperature, y=rate.output)) + geom_line()
ggplot(posterior_draws, aes(x=Temperature, y=.linpred, color=.linpred)) +
geom_point(alpha=0.05, size=1) +
geom_smooth(data=as.data.frame(spline_fit$data), aes(x=x, y=y), color="red",
size=3) + ylab("rate.output")+
geom_vline(xintercept = my.seg$psi[1,2]) +
geom_vline(xintercept = my.seg$psi[2,2])+
geom_vline(xintercept = my.seg$psi[3,2])+
geom_vline(xintercept = my.seg$psi[1,1], color="blue") +
geom_vline(xintercept = my.seg$psi[2,1], color="blue") +
geom_vline(xintercept = my.seg$psi[3,1], color="blue") +
geom_line(data=my.model, aes(x=temperature, y=rate.output), color="black", size=2) +
theme_classic()## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'